/************************************************************************************
*   Project:        Angina paper
*	Goal:			Summary Statistics and Main Regressions

*   Input Data:     Angina_mainregs.dta
*                   Angina_vcdata

*   Last Modified: 	2020-08-08

**************************************************************************************/

// Summary Statistics

******************************************************************
*             Clinic Characteristics (Table 2)                   *
******************************************************************

use "$data/working/Angina_vcdata.dta", clear

	table1, vars(fd_inc contn %9.2f \ fd_clin contn %9.2f \ fd_nclin contn %9.2f \ fd_pat contn %9.2f \ ///
		fd_typemed contn %9.2f \ fd_typewmed contn %9.2f \ fd_typecmed contn %9.2f \ fd_valinst contn %9.2f \ ///
		fd_popu_ph contn %9.2f \ fd_instele bine \ fd_instste bine \ fd_instsph bine \ fd_instthe bin) cformat(%12.0f) ///
		saving("$output/summary_clinic1.xls",replace) 

	global VC1 fd_inc fd_clin fd_nclin fd_pat fd_typemed fd_typewmed fd_typecmed fd_valinst fd_popu_ph  // continuous variables
	global VC2 fd_instele fd_instste fd_instsph fd_instthe                                              // binary variables

	d $VC1 $VC2	
	ci means $VC1
	ci proportions $VC2, wilson

	eststo clear
		
		eststo: estpost ci $VC1
		eststo: estpost ci $VC2, binomial wilson                    
		esttab, cells("b(fmt(2)) se(fmt(2)) lb(fmt(2)) ub(fmt(2))") label

		esttab using "$output/summary_clinic2.csv", cells("b(fmt(2)) se(fmt(2)) lb(fmt(2)) ub(fmt(2))") label replace


	* Clinic Characteristics by Provinces
	gen PID=substr(vcid,1,1)
	destring PID, replace
	label var PID "Province ID"

	table1, by(PID) vars(fd_inc contn %9.2f \ fd_clin contn %9.2f \ fd_nclin contn %9.2f \ fd_pat contn %9.2f \ ///
		fd_typemed contn %9.2f \ fd_typewmed contn %9.2f \ fd_typecmed contn %9.2f \ fd_valinst contn %9.2f \ ///
		fd_popu_ph contn %9.2f \ fd_instele bin \ fd_instste bin \ fd_instsph bin \ fd_instthe bin) cformat(%12.0f) ///
		saving("$output/summary_clinic3.xls",replace) 	


******************************************************************
*         Village Clinician Characteristics (Table 1)            *
******************************************************************

use "$data/working/Angina_mainregs.dta", clear

	table1, vars(vd_age contn %9.2f \ vd_sex bin \ vd_exp contn %9.2f \ vd_hied cat \ ///
		vd_himed cat \ vd_typmed cat \ vd_qual cat \ vd_bsal_binary bin \ vd_tinc contn %9.2f \ vd_hours contn %9.2f) cformat(%12.0f) ///
		saving("$output/summary_doctor1.xls",replace)

		foreach v in vd_hied vd_himed vd_typmed vd_qual {
			tab `v', gen(`v'_) 
			}

	global VD1 "vd_age vd_exp vd_tinc vd_hours"                                                                  // continuous
	global VD2 "vd_sex vd_hied_1-vd_hied_3 vd_himed_1-vd_himed_3 vd_typmed_1-vd_typmed_5 vd_qual_1-vd_qual_4 vd_bsal_binary"   // binary

	d $VD1 $VD2		
	ci means $VD1 
	ci proportions $VD2, wilson

	eststo clear
		
		eststo: estpost ci $VD1
		eststo: estpost ci $VD2, binomial wilson                    
		esttab, cells("b(fmt(2)) se(fmt(2)) lb(fmt(2)) ub(fmt(2))") label

		esttab using "$output/summary_doctor2.csv", cells("b(fmt(2)) se(fmt(2)) lb(fmt(2)) ub(fmt(2))") label replace


	* Doctor Characteristics by qualification (Appendix)
	tab vd_qual
	gen vd_qual_hi = (vd_qual>2)
	replace vd_qual_hi=. if vd_qual==.
	label define vd_qual_hi 0 "[0] Village doctor qualification or lower" 1 "[1] Higher than village doctor qualification"
	label val vd_qual_hi vd_qual_hi
	label var vd_qual_hi "Higher than village doctor qualification (0/1)"
	tab vd_qual_hi

	table1, by(vd_qual_hi) vars(vd_age contn %9.2f \ vd_sex bin \ vd_exp contn %9.2f \ vd_hied cat \ ///
		vd_himed cat \ vd_typmed cat \ vd_bsal_binary bin \ vd_tinc contn %9.2f \ vd_hours contn %9.2f) cformat(%12.0f) ///
		saving("$output/summary_doctor3.xls",replace)	


	* Doctor Characteristics by Provinces (Appendix)
	gen PID=substr(vcid,1,1)
	destring PID, replace
	label var PID "Province ID"

	table1, by(PID) vars(vd_age contn %9.2f \ vd_sex bin \ vd_exp contn %9.2f \ vd_hied cat \ ///
		vd_himed cat \ vd_typmed cat \ vd_qual cat \ vd_bsal_binary bin \ vd_tinc contn %9.2f \ vd_hours contn %9.2f) cformat(%12.0f) ///
		saving("$output/summary_doctor4.xls",replace)


******************************************************************
*            Recommended Questions (Table A1)                    *
******************************************************************
	
	table1, vars(rec_q1 bine \ rec_q2 bine \ rec_q3 bine \ rec_q4 bine \ rec_q5 bine \ ///
		rec_q6 bine \ rec_q7 bine \ rec_q8 bine \ rec_q9 bine \ rec_q10 bine \ rec_q11 bine \ ///
		rec_q13 bine \ rec_q14 bine \ rec_q15 bine \ rec_q16 bine ) cformat(%12.0f) ///
		saving("$output/summary_recques1.xls",replace)                              // not include rec_q12=0


	global rec_q rec_q1 rec_q2 rec_q3 rec_q4 rec_q5 rec_q6 rec_q7 rec_q8 rec_q9 ///
		rec_q10 rec_q11 rec_q12 rec_q13  rec_q14 rec_q15 rec_q16      
	global ess_q ess_q1 ess_q2 ess_q3 ess_q4 ess_q5 ess_q6 ess_q7

	d $rec_q $ess_q
	ci proportions $rec_q $ess_q, wilson

	eststo clear

		eststo: estpost ci $rec_q $ess_q, binomial wilson                    		
		esttab, cells("b(fmt(2)) se(fmt(2)) lb(fmt(2)) ub(fmt(2))") label

		esttab using "$output/summary_recques2.csv", cells("b(fmt(2)) se(fmt(2)) lb(fmt(2)) ub(fmt(2))") nomtitle nonumber label replace


******************************************************************
*             Recommended Exams (Table A2)                       *
******************************************************************

	table1, vars(rec_e1 bine \ rec_e2 bine \ rec_e3 bine \ rec_e4 bine \ rec_e5 bine \ rec_e6 bine) cformat(%12.0f) ///
		saving("$output/summary_recexam1.xls",replace)

	global rec_e rec_e1 rec_e2 rec_e3 rec_e4 rec_e5 rec_e6 
	global ess_e ess_e1 

	d $rec_e $ess_e

	ci proportions $rec_e $ess_e, wilson

	eststo clear

		eststo: estpost ci $rec_e $ess_e, binomial wilson                    		
		esttab, cells("b(fmt(2)) se(fmt(2)) lb(fmt(2)) ub(fmt(2))") label

		esttab using "$output/summary_recexam2.csv", cells("b(fmt(2)) se(fmt(2)) lb(fmt(2)) ub(fmt(2))") nomtitle nonumber label replace


******************************************************************
*    Table 3  Competence in Process and Diagnosis (Table 3)      *
******************************************************************

	table1, vars(asked_q contn %9.2f \ rec_q_num contn %9.2f \ rec_q_percent contn %9.2f \ ess_q_num contn %9.2f \ ///
		ess_q_percent contn %9.2f \ performed_e contn %9.2f \ rec_e_num contn %9.2f \ rec_e_percent contn %9.2f \ ///
		ess_e_num contn %9.2f \ ess_e_percent contn %9.2f \ rec_qande_num contn %9.2f \ rec_qande_percent contn %9.2f \ ///
		ess_qande_num contn %9.2f \ ess_qande_percent contn %9.2f) ///
		saving("$output/summary_process1.xls",replace)

	table1, vars(v_gavediagnosis bine \ v_diagnosis cat) cformat(%12.0f) saving("$output/summary_diagnosis1.xls",replace)


	global process asked_q rec_q_num rec_q_percent ess_q_num ess_q_percent ///
		performed_e rec_e_num rec_e_percent ess_e_num ess_e_percent ///
		rec_qande_num rec_qande_percent ess_qande_num ess_qande_percent 

	global diagnosis v_gavediagnosis v_corrdiag v_partdiag v_wrongdiag

	d $process $diagnosis	
	ci means $process
	ci proportions $diagnosis, wilson

	eststo clear
		
		eststo: estpost ci $process
		eststo: estpost ci $diagnosis, binomial wilson                    
		esttab, cells("b(fmt(2)) se(fmt(2)) lb(fmt(2)) ub(fmt(2))") label

		esttab using "$output/summary_process2_diagnosis2.csv", cells("b(fmt(2)) se(fmt(2)) lb(fmt(2)) ub(fmt(2))") label replace


******************************************************************
*               Management of Angina (Table 4)                   *
******************************************************************

use "$data/working/Angina_mainregs.dta", clear

	keep vcdocid q6_* q7_* q10_* q11_*
	drop if q7_correcttreat==. & q11_correcttreat==.   // drop 2 doctors who did not treat

		foreach var of varlist q6_* q7* {
			local newname=substr("`var'", 4, .)
	   		rename `var' `newname'0
		}

		foreach var of varlist q10_* q11* {
			local newname=substr("`var'", 5, .)
	   		rename `var' `newname'1
		}
		rename v_gavetreat1 gavetreat1

		reshape long correcttreat correcttreat_1 correcttreat_2 partcorrecttreat partcorrecttreat_1 partcorrecttreat_2 ///
			partcorrecttreat_3 partcorrecttreat_4 partcorrecttreat_5 incorrecttreat incorrecttreat_1 incorrecttreat_2 incorrecttreat_3 ///
			gavereferral gavetreat numofdrugs numofwestern numofchn numofherb correctmed unnecessary unnecessary_west ///
			unnecessary_chn unnecessary_herb harmmed antibio hormone analgesic, i(vcdocid) j(truedia)

		label var truedia "=1 if true diagnosis given"
		label var correcttreat "=1 if correct treatment"
		label var correcttreat_1 "=1 if correct treatment 1: referral & no drugs"
		label var correcttreat_2 "=1 if correct treatment 2: referral & correct drugs & no harmful/unnecessary"
		label var partcorrecttreat "=1 if partially correct treatment"
		label var partcorrecttreat_1 "=1 if partially correct 1: referral & correct drugs & harmful/unnecessary"
		label var partcorrecttreat_2 "=1 if partially correct 2: referral & harmful drugs"
		label var partcorrecttreat_3 "=1 if partially correct 3: referral & unnecessary drugs"
		label var partcorrecttreat_4 "=1 if partially correct 4: no referral & correct drugs & no harmful/unnecessary"		
		label var partcorrecttreat_5 "=1 if partially correct 5: no referral & correct drugs & harmful/unnecessary"		
		label var incorrecttreat "=1 if incorrect treatment"
		label var incorrecttreat_1 "=1 if incorrect treatment 1: no referral & no drugs"
		label var incorrecttreat_2 "=1 if incorrect treatment 2: no referral & harmful drugs"
		label var incorrecttreat_3 "=1 if incorrect treatment 3: no referral & unnecessary drugs"
		label var gavereferral "=1 if patient was referred"
		label var gavetreat "=1 if gave any medication"
		label var numofdrugs "Number of drugs prescribed"
		label var numofwestern "Number of western drugs"
		label var numofchn "Number of Chinese patent drugs"
		label var numofherb "Number of herbal drugs"
		label var correctmed "=1 if correct medication"
		label var unnecessary "=1 if unnecessary medication"
		label var unnecessary_west "=1 if unnecessary western medication"
		label var unnecessary_chn "=1 if unnecessary Chinese patent medication"
		label var unnecessary_herb "=1 if unnecessary herbal medication"
		label var harmmed "=1 if harmful medication"
		label var antibio "=1 if any antibiotics prescribed"
		label var hormone "=1 if any hormone prescribed"
		label var analgesic "=1 if any analgesic prescribed"

		* Province ID	
		gen PID=substr(vcdocid,1,1)
		destring PID, replace
		label var PID "Province ID"

	table1, by(truedia) vars(correcttreat bine \ correcttreat_1 bine \ correcttreat_2 bine \ partcorrecttreat bine \ ///
		partcorrecttreat_2 bine \ partcorrecttreat_3 bine \ partcorrecttreat_4 bine \ ///
		incorrecttreat bine \ incorrecttreat_1 bine \ incorrecttreat_2 bine \ incorrecttreat_3 bine \ gavereferral bine \ ///
		gavetreat bine \ numofdrugs contn %9.2f \ numofwestern contn %9.2f \ numofchn contn %9.2f \ numofherb contn %9.2f \ ///
		correctmed bine \ unnecessary bine \ unnecessary_west bine \ unnecessary_chn bine \  ///
		harmmed bine \ antibio bine \ hormone bine \ analgesic bine) cformat(%12.0f) ///
		saving("$output/summary_management1.xls",replace)       // not include unnecessary_herb partcorrecttreat_1 partcorrecttreat_5 (all 0)

	* Means & CIs
	global management correcttreat correcttreat_1 correcttreat_2 partcorrecttreat partcorrecttreat_1 partcorrecttreat_2 ///
		partcorrecttreat_3 partcorrecttreat_4 partcorrecttreat_5 incorrecttreat incorrecttreat_1 incorrecttreat_2 incorrecttreat_3
	global component1 gavereferral gavetreat correctmed unnecessary unnecessary_west unnecessary_chn unnecessary_herb ///
		harmmed antibio hormone analgesic
	global component2 numofdrugs numofwestern numofchn numofherb    // continous 

	d $management $component1 $component2	

	eststo clear

		eststo: estpost ci $management if truedia==0, binomial wilson                    		
		eststo: estpost ci $component1 if truedia==0, binomial wilson                    
		eststo: estpost ci $component2 if truedia==0

		eststo: estpost ci $management if truedia==1, binomial wilson                    		
		eststo: estpost ci $component1 if truedia==1, binomial wilson                    
		eststo: estpost ci $component2 if truedia==1

		esttab, cells("b(fmt(2)) se(fmt(2)) lb(fmt(2)) ub(fmt(2))") label

		esttab using "$output/summary_management2.csv", cells("b(fmt(2)) se(fmt(2)) lb(fmt(2)) ub(fmt(2))") label replace


	global management correcttreat correcttreat_1 correcttreat_2 partcorrecttreat partcorrecttreat_2  ///
		partcorrecttreat_3 partcorrecttreat_4 incorrecttreat incorrecttreat_1 incorrecttreat_2 incorrecttreat_3
	global component1 gavereferral gavetreat correctmed unnecessary unnecessary_west unnecessary_chn  ///
		harmmed antibio hormone analgesic
	global component2 numofdrugs numofwestern numofchn numofherb   
 
	foreach v in $management $component1 {                   // not include unnecessary_herb partcorrecttreat_1 partcorrecttreat_5 (all 0)             
		tab `v'
		prtest `v', by(truedia)
		}	

	foreach v in $component2 {		
		sum `v'
		ttest `v', by(truedia)
		}		


// Main Regressions

use "$data/working/Angina_mainregs.dta", clear 

** Data Preparation

	gen PID=substr(vcdocid,1,1)
	destring PID, replace
	label var PID "Province ID"

		* Recode type of medical education
		tab vd_typmed 
		recode vd_typmed 1 5=1 2=2 3=3 4=4, gen(vd_typmed2)            // recode pulic health as no medical education
		label define vd_typmed2 1 "[1] No Medical Education" 2 "[2] Chinese" 3 "[3] Western" 4 "[4] Chinese & Western"
		label val vd_typmed2 vd_typmed2
		tab vd_typmed2

		* rescale 
		sum vd_tinc
		gen vd_tinck=vd_tinc/1000
		label var vd_tinck "Total annual income of village doctors (in thousands)"
		sum vd_tinck

			gen vd_tinc_ln = ln(vd_tinc)
			label var vd_tinc_ln "ln(total annual income of village doctor)"
			sum vd_tinc_ln

		sum fd_inc
		gen fd_inck=fd_inc/1000
		label var fd_inck "Annual net income of village clinics (in thousands)"
		sum fd_inck

			gen fd_inc_ln = ln(fd_inc)
			label var fd_inc_ln "ln(Annual net income of village clinics)"

		sum fd_pat
		gen fd_patln=ln(fd_pat)
		label var fd_patln "ln(annual patient volume)"
		sum fd_patln

		sum fd_typemed
		gen fd_typemedln=ln(fd_typemed)
		label var fd_typemedln "ln(varieties of medicine in inventory)"
		sum fd_typemedln

		sum fd_valinst
		gen fd_valinstk=fd_valinst/1000
		label var fd_valinstk "Total value of medical instruments (in thousands)"
		sum fd_valinstk

			gen fd_valinst_ln = ln(fd_valinst)
			label var fd_valinst_ln "ln(total value of medical instruments)"

		sum fd_popu_ph
		gen fd_popu_phk=fd_popu_ph/1000
		label var fd_popu_phk "pop for public health services (in 1000)"
		sum fd_popu_phk

		label var vd_sex "male village doctor (0/1)"


	*Define binary variables
	g vd_ed_high = (vd_hied>2)
		replace vd_ed_high=. if vd_hied==. 
		label var vd_ed_high "higher education (0/1)"

	g nomeded = (vd_typmed2==1)
		replace nomeded =. if vd_typmed2==.

	g vd_qual_hi = (vd_qual>2)
		replace vd_qual_hi=. if vd_qual==.
		label var vd_qual_hi "higher qualification (0/1)"

	su vd_exp, det
	gen vd_exp_hi = (vd_exp>`r(p50)')
		replace vd_exp_hi=. if vd_exp==.
		label var vd_exp_hi "longer practicing years (0/1)"

	global VD "vd_sex vd_exp_hi vd_ed_high vd_qual_hi vd_tinc_ln" 
	global VC "fd_clin fd_patln fd_typemedln fd_valinst_ln fd_popu_phk" 

	sum $VD $VC

******************************************************************
*        Regression Analysis I -- Diagnostics                    *
******************************************************************

* Generate correct or partially correct
	gen v_corrpartdiag=(v_corrdiag==1|v_partdiag==1)
	label var v_corrpartdiag "correct and partially correct diagnosis"

	global Y1 "rec_qande_percent"        // prop of recommended questions and exams
	global Y2 "ess_qande_percent"		// prop essential items
 	global Y3 "v_corrpartdiag"           // correct or partially correct diagnosis

 	d $Y1 $Y2 $Y3
	sum $Y1 $Y2 $Y3                        

	* Non-linear
	eststo clear

		* Regression on the probability of a correct or partially correct diagnosis
		logit $Y3 $VD $VC i.PID, vce(cluster vcid)
		eststo, title(Diag): margins, dydx(*) post		

		* Regression on the proportion of recommended questions and examinations completed
		glm $Y1 $VD $VC i.PID, family(binomial) link(logit) vce(cluster vcid)
		eststo, title(Rec): margins, dydx(*) post

		* Regression on the proportion of essential questions and examinations completed
		glm $Y2 $VD $VC i.PID, family(binomial) link(logit) vce(cluster vcid)  
		eststo, title(Ess): margins, dydx(*) post

		esttab, label b(2) ci(2) star(* 0.10 ** 0.05 *** 0.01) replace nogaps mtitles

		esttab using "$output/DiagnosticReg.csv", b(%9.2fc) ci(%9.2fc) ///
			starlevels( * 0.1 ** 0.05 *** 0.01) ar2(2) label replace nobase nogap mti compress  ///
			title("Correlates of Diagnostic Quality")

******************************************************************
*                      Regressions -- Treatment                  *
******************************************************************

*	drop if q7_correcttreat==. & q11_correcttreat==.   // drop 2 doctors who did not treat

	* generate correct/partially treatment
	gen q7_corrparttreat=(q7_correcttreat==1|q7_partcorrecttreat==1)
		replace q7_corrparttreat=. if q7_correcttreat==. & q11_correcttreat==.
		label var q7_corrparttreat "correct and partially correct treatment"	

	*** Regressions for managment under own diagnosis

		gen q6_unneharm=(q6_unnecessary==1|q6_harmmed==1) if q6_gavetreat==1
		label var q6_unneharm "unnecessary and harmful drugs"

	global Y11 "q7_corrparttreat"        // correct and partially correct treatment
	global Y12 "q7_gavereferral"         // referral
	global Y13 "q6_gavetreat"            // any medication
	global Y14 "q6_numofdrugs"           // number of drugs
	global Y15 "q6_correctmed"           // correct medication
	global Y16 "q6_unneharm"             // unnecessary and harmful drugs

	d $Y11 $Y12 $Y13 $Y14 $Y15 $Y16
	sum $Y11 $Y12 $Y13 $Y14 $Y15 $Y16

	* Treatment - Own diagnosis
	eststo clear

		* Correct and partially correct treatment
		logit $Y11 $VD $VC i.PID, vce(cluster vcid)
		margins, dydx(*) post

			// perfect prediction for one prov, observations drop --> Firth logit

			capture prog drop firthmarg
			prog firthmarg, eclass
				firthlogit  $Y11 $VD $VC i.PID
				margins, dydx(*) expression(invlogit(predict(xb))) post
			end 

			eststo, title(cpctreat): bootstrap _b, cluster(vcid) reps(400): firthmarg
			

		* Referral	
		logit $Y12 $VD $VC i.PID, vce(cluster vcid) 
		margins, dydx(*) post	

			capture prog drop firthmarg
			prog firthmarg, eclass
				firthlogit  $Y12 $VD $VC i.PID
				margins, dydx(*) expression(invlogit(predict(xb))) post
			end 

			eststo, title(referral): bootstrap _b, cluster(vcid) reps(400): firthmarg


		* Any medication
		logit $Y13 $VD $VC i.PID, vce(cluster vcid) 
		eststo, title(anymed): margins, dydx(*) post

		* Number of drugs
		reg $Y14 $VD $VC i.PID, vce(cluster vcid) 
		eststo, title(numdrugs): margins, dydx(*) post

		* Harmful medication
		logit $Y16 $VD $VC i.PID, vce(cluster vcid) 
		eststo, title(unnecharmmed): margins, dydx(*) post

	esttab, label b(2) ci(2) star(* 0.10 ** 0.05 *** 0.01) replace nogaps mtitles

	esttab using "$output/TreatReg.csv", b(%9.2fc) ci(%9.2fc) ///
			starlevels( * 0.1 ** 0.05 *** 0.01) ar2(2) label replace nobase nogap mti compress  ///
			title("Correlates of Correct Treatment")	


	* Treatment - True diagnosis
		gen q11_corrparttreat=(q11_correcttreat==1|q11_partcorrecttreat==1)
		label var q11_corrparttreat "correct and partially correct treatment"

		gen q10_unneharm=(q10_unnecessary==1|q10_harmmed==1) if q10_v_gavetreat==1
		label var q10_unneharm "unnecessary and harmful drugs"

	global Y21 "q11_corrparttreat"          // correct and partially correct treatment
	global Y22 "q11_gavereferral"           // referral
	global Y23 "q10_v_gavetreat"            // any medication
	global Y24 "q10_numofdrugs"             // number of drugs
	global Y25 "q10_correctmed"             // correct medication
	global Y26 "q10_unneharm"               // unnecessary and harmful drugs

	d $Y21 $Y22 $Y23 $Y24 $Y25 $Y26
	sum $Y21 $Y22 $Y23 $Y24 $Y25 $Y26

	* Non-linear
	eststo clear

		* Correct and partially correct treatment
		logit $Y21 $VD $VC i.PID, vce(cluster vcid)
		eststo, title(treat2): margins, dydx(*) post

		* Referral	
		logit $Y22 $VD $VC i.PID, vce(cluster vcid) 
		margins, dydx(*) 	

			capture prog drop firthmarg
			prog firthmarg, eclass
				firthlogit  $Y22 $VD $VC i.PID
				margins, dydx(*) expression(invlogit(predict(xb))) post
			end 

			eststo, title(referral2): bootstrap _b, cluster(vcid) reps(400): firthmarg

		* Any medication
		logit $Y23 $VD $VC i.PID, vce(cluster vcid) 
		eststo, title(anymed2): margins, dydx(*) post

		* Number of drugs
		reg $Y24 $VD $VC i.PID, vce(cluster vcid) 
		eststo, title(numdrugs2): margins, dydx(*) post

		* Harmful medication
		logit $Y26 $VD $VC i.PID, vce(cluster vcid) 
		eststo, title(unnecharmmed2): margins, dydx(*) post

	esttab, label b(2) ci(2) star(* 0.10 ** 0.05 *** 0.01) replace nogaps mtitles

	esttab using "$output/TreatReg.csv", b(%9.2fc) ci(%9.2fc) ///
			starlevels( * 0.1 ** 0.05 *** 0.01) ar2(2) label append nobase nogap mti compress  ///
			title("Correlates of Correct Treatment under 'True' Diagnosis")	


******************************************************************
*            Diagnosis, TREATMENT AND CHECKLIST COMPLETION       *
******************************************************************

*OWN DIAGNOSIS
g corrmed = (q6_correctmed == 1) // unconditional correct medication

	eststo clear

			logit v_corrpartdiag rec_qande_percent $VD $VC i.PID, vce(cluster vcid)
			eststo: margins, dydx(*) post

			logit q7_corrparttreat rec_qande_percent $VD $VC i.PID, vce(cluster vcid)
			margins, dydx(*) 

					capture prog drop firthmarg
					prog firthmarg, eclass
						firthlogit  q7_corrparttreat rec_qande_percent $VD $VC i.PID
						margins, dydx(*) expression(invlogit(predict(xb))) post
					end 

					eststo: bootstrap _b, cluster(vcid) reps(400): firthmarg
			
			logit q7_gavereferral rec_qande_percent $VD $VC i.PID, vce(cluster vcid)
			margins, dydx(*) 

					capture prog drop firthmarg
					prog firthmarg, eclass
						firthlogit  q7_gavereferral rec_qande_percent $VD $VC i.PID
						margins, dydx(*) expression(invlogit(predict(xb))) post
					end 

					eststo: bootstrap _b, cluster(vcid) reps(400): firthmarg

			logit corrmed rec_qande_percent $VD $VC i.PID, vce(cluster vcid)
			eststo: margins, dydx(*) post

			logit q6_unneharm rec_qande_percent $VD $VC i.PID, vce(cluster vcid)
			eststo: margins, dydx(*) post

		esttab, label b(2) ci(2) star(* 0.10 ** 0.05 *** 0.01) replace nogaps mtitles

		esttab using "$output/Treat-CL.csv", b(%9.2fc) ci(%9.2fc) ///
				starlevels( * 0.1 ** 0.05 *** 0.01) ar2(2) label replace nobase nogap mti compress  ///
				title("Correlation Between Correct Treatment and Checklist Completion") ///
				keep(rec_qande_percent)

*TRUE DIAGNOSIS
g corrmed_true = (q10_correctmed==1) 


	eststo clear

			logit v_corrpartdiag rec_qande_percent $VD $VC i.PID, vce(cluster vcid) // placeholder
			eststo: margins, dydx(*) post

			logit q11_corrparttreat rec_qande_percent $VD $VC i.PID, vce(cluster vcid)
			eststo: margins, dydx(*) post

			logit q11_gavereferral rec_qande_percent $VD $VC i.PID, vce(cluster vcid)
			margins, dydx(*) 
					
					capture prog drop firthmarg
					prog firthmarg, eclass
						firthlogit  q11_gavereferral rec_qande_percent $VD $VC i.PID
						margins, dydx(*) expression(invlogit(predict(xb))) post
					end 

					eststo: bootstrap _b, cluster(vcid) reps(400): firthmarg

			logit corrmed_true rec_qande_percent $VD $VC i.PID, vce(cluster vcid)
			eststo: margins, dydx(*) post

			logit q10_unneharm rec_qande_percent $VD $VC i.PID, vce(cluster vcid)
			eststo: margins, dydx(*) post

		esttab, label b(2) ci(2) star(* 0.10 ** 0.05 *** 0.01) replace nogaps mtitles

		esttab using "$output/Treat-CL.csv", b(%9.2fc) ci(%9.2fc) ///
			starlevels( * 0.1 ** 0.05 *** 0.01) ar2(2) label append nobase nogap mti compress  ///
			title("Correlation Between Correct Treatment and Checklist Completion - True Diagnosis") ///
			keep(rec_qande_percent)

exit
